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Cell biologists increasingly rely on computer-aided 
image analysis, allowing them to collect precise, 
unbiased quantitative results. However, despite 
great progress in image processing and computer vision, 
current computational approaches fail to address many 
key aspects of cell behavior, including the cell protrusions 
that guide cell migration and drive morphogenesis. We 
developed the open source MATLAB application CellGeo, 
a user-friendly computational platform to allow simultane- 
ous, automated tracking and analysis of dynamic changes 
in cell shape, including protrusions ranging from filopodia 



to lamellipodia. Our method maps an arbitrary cell shape 
onto a tree graph that, unlike traditional skeletonization 
algorithms, preserves complex boundary features. CellGeo 
allows rigorous but flexible definition and accurate auto- 
mated detection and tracking of geometric features of in- 
terest. We demonstrate CellGeo's utility by deriving new 
insights into (a) the roles of Diaphanous, Enabled, and 
Capping protein in regulating filopodia and lamellipodia 
dynamics in Drosophila melanogaster cells and (b) the 
dynamic properties of growth cones in catecholaminergic a- 
differentiated neuroblastoma cells. 



Introduction 

Cell protrusions are an essential driver of dynamic cell shape 
changes and motility during development and disease. Morpho- 
genic processes from gastrulation to organogenesis require co- 
ordinated protrusive behavior to shape tissues and organs. Cell 
protrusions are also essential for cell migration during wound 
healing, and cancer cells use protrusions to migrate from pri- 
mary tumors during metastasis. Cells use both lamellipodia and 
filopodia to generate shape changes and drive motility, making 
it imperative to understand how the dynamics of both structures 
are regulated. 

Recent advances in live-cell imaging, including new mi- 
croscope designs and novel molecular probes, allowed biolo- 
gists to visualize cellular behavior with extraordinary precision 
and detail. However, to fully take advantage of these advances 
requires novel computational methods for image processing 
and analysis (Meijering et al., 2004; Costantino et al., 2008; 
Fanti et al., 201 1). Here, we present the computational platform 
CellGeo, a MATLAB application to identify, track, and charac- 
terize dynamic cell shape changes (Fig. 1 A). The key step in 

Correspondence to Timothy C. Elston: timothy_elston@med.unc.edu; or Mark 
Peifer: peifer@unc.edu 

Abbreviations use in this paper: CAD, catecholaminergic a differentiated; CP, 
Capping protein; DAD, Dia autoinhibitory domain; Dia, Diaphanous; Ena, En- 
abled; GUI, graphical user interface; LPA, lysophosphatidic acid; MAT, medial 
axis transform; mCh, monomeric Cherry. 



CellGeo is the representation of any arbitrary cell shape as a 
tree graph (Fig. 1, C-F; and Video 1). This conversion facili- 
tates precise definitions of shape features, such as filopodia and 
lamellipodia, and quantitative analyses of their dynamics. Cell- 
Geo is a fully automated system with a graphical user interface 
(GUI) for easy adjustment of parameters for flexible and accu- 
rate protrusion and cell body detection and analysis of any cell 
type (Fig. 1 A). CellGeo has an intuitive/self-explanatory design 
that allows two modes of operation: (1) an interactive exploratory 
mode, in which users can see how changes in parameters affect 
the analysis and adjust them accordingly; and (2) an unsuper- 
vised production mode, in which users simply import data, click 
a button, and save results using default or previously set param- 
eter values. 

To demonstrate CellGeo's utility, we use the platform to 
analyze the effects of the actin regulators Diaphanous (Dia), En- 
abled (Ena), and Capping protein (CP) on both filopodial and la- 
mellipodial dynamics in Drosophila melanogaster cells (Fig. 1 B). 
Ena and Dia both localize to filopodia and lamellipodia, and 
overexpressing either drives both types of protrusions. However, 
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Figure 1 . CellGeo platform architecture and qualitative interpretation of the MAT. (A) CellGeo package pipeline for defining, detecting and tracking both 
thin or broad cellular protrusions or growth cones. (B) Drosophila Dl 6C3 cell expressing GFP-actin with four kymographs labeled 1-4 show high variability 
of protrusiveness within a single cell, making analysis biased by placement. Bar, 5 pm. (C) A cell and its boundary (yellow). Bar, 10 urn. (D] Distance 
function F(x, y), which measures the minimum distance from each point (x, y) inside the cell to the cell boundary. (E) Ridges (surface discontinuities) of the 
minimum distance function in D. (F) Medial axis as a top view of the ridges in E. 



only Ena is required for filopodia number and length, demon- 
strating that CellGeo can identify distinct roles of similar actin 
regulators in controlling the complex composition of cell pro- 
trusions (Videos 2-6). We also find Ena and Dia can act inde- 
pendently of one another in the formation of filopodia and broad 
protrusions and that CP is required to limit Ena activity, likely 



by limiting availability of barbed ends. To further demonstrate 
CellGeo's versatility, we use it to study neuronal growth cone dy- 
namics and the role of the GTPase RhoA (Etienne-Manneville 
and Hall, 2002; Jaffe and Hall, 2005) in driving this behavior. 
Our analysis revealed in a quantitative way the spatiotemporal 
distribution of RhoA activity in growth cones and cell bodies 
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during growth cone protrusion and retraction. Our analysis also 
revealed unexpected correlations between geometric character- 
istics of growth cones and the delay in onset of growth cone re- 
traction after treatment with lysophosphatidic acid (LPA). 

Results 

Defining cell edge features 

The complex structure and dynamics of cellular protrusions 
reflects their functional diversity, making designing automated 
algorithms for quantifying their behavior challenging. We de- 
veloped a novel approach for identifying and tracking protru- 
sions that functions reliably even in the presence of irregular 
geometries and complex dynamics. 

Our method begins by identifying the cell boundary in 
each frame of a time-lapse video. Although other commonly 
used software packages like ImageJ or MetaMorph provide seg- 
mentation routines, we included a simple segmentation module, 
MovThresh, to incorporate additional features for dynamic shape 
analysis and make CellGeo a stand-alone platform. MovThresh 
automatically thresholds each frame of a video using the Otsu 
algorithm or one of the other built-in segmentation routines, 
overlays the detected boundary with the original image, and dis- 
plays suggested threshold values as a function of time (Fig. SI). 
If the algorithm accurately captured the cell edge in each frame, 
no additional steps are required. Often, however, the threshold 
values determined by the algorithm vary from frame to frame. 
Therefore, depending on the type of analysis to be performed, 
users might opt to use a smooth curve through the suggested 
threshold values or a constant (e.g., mean) value to rethreshold 
all frames consistently. In special cases, when image brightness 
significantly changes during the experiment (e.g., induction of a 
fluorescent protein) users can interactively set the threshold in a 
way appropriate for the current analysis. The main purpose of 
these steps is to avoid or minimize boundary artifacts, such as 
unresolved or merged filopodia, that result from over- or under- 
thresholding images (Fig. 2, A-C). 

The output of MovThresh is a closed polygon made from 
lines connecting the centers of neighboring boundary pixels 
(Fig. S2). In addition to being computationally efficient and 
simple, this representation preserves the pixelated nature of the 
raw data and does not require any image preprocessing. Impor- 
tantly, the next CellGeo module, BisectoGraph, is based on a 
grid-free approach so that no association with pixel centers is 
necessary: BisectoGraph can process any simple closed poly- 
gon specified by (x, y) coordinates. 

BisectoGraph is the key step in analyzing cell edge dy- 
namics: it creates a representation of the entire cell that makes 
quantifying protrusions and cell body dynamics precise and un- 
biased. BisectoGraph relies on a Voronoi transformation that 
converts the cell boundary found by MovThresh into a tree 
graph. Voronoi diagrams were originally introduced as a way to 
partition space into regions that uniquely enclose individual 
points of an arbitrary discrete set (Aurenhammer, 1991). How- 
ever, the Voronoi transformation also can be applied to a dis- 
crete set of line segments or, as a limiting case, to a closed 
polygon. The Voronoi diagram inside a polygon is a union of the 



medial axis and a set of rays originating at concave boundary 
points (Materials and methods; Fig. S4; Mayya and Rijan, 1995). 

The concept of the medial axis transform (MAT; or sym- 
metric axis or skeleton) was introduced by Blum (1967) as a 
method for cell shape analysis. The MAT of a plane shape is de- 
fined as the set formed from the union of the centers of all maxi- 
mal inscribed circles (circles contained within the shape that share 
at least two boundary points with it; Chin et al., 1999). For ex- 
ample, the MAT of a rectangle consists of five branches (Fig. 3 A). 
Since its introduction, the MAT found applications in many 
areas, including object recognition, graphics, 3D modeling, 
computer-aided manufacturing, and characterizing structures 
of proteins or organs (Morse et al., 1994; Styner et al., 2003; 
Abeysinghe et al., 2008). 

A common theme in all these applications is the genera- 
tion of a reduced representation of the shape of interest that 
captures its large-scale features. Therefore, algorithms for con- 
structing MATs are typically accompanied by pruning algorithms 
that decide which medial axis branches can safely be removed 
(Tek and Kimia, 2003; Aichholzer et al., 2009). A major short- 
coming of this is that small perturbations to shape can generate 
disproportionately large perturbations to the MAT, necessitating 
postprocessing to reduce the effects of noise. In contrast, we 
demonstrate that for the analysis of protrusion dynamics, it is 
useful to keep the complete Voronoi tree graph with branches 
reaching every boundary point and define boundary measures 
that are insensitive to noise (Materials and methods). Fig. 2 
illustrates the fact that although small (less than 1 pixel) bound- 
ary perturbations produce significant differences in the tree 
graph (Fig. 2 D, white and black lines), the detection of protru- 
sion tips using these graphs by means of the boundary profile is 
not significantly affected. Indeed, Fig. 2 (E-G) shows 100 over- 
laid boundaries and the detected tips for each boundary (Fig. 2, 
E-G, red dots), whereas Fig. 2 (H-I) shows the resulting 
boundary profiles and the local maxima corresponding to the 
tips (Fig. 2, H-I, also marked by red dots). 

To illustrate the effect of small perturbations, consider a 
digital image of an inclined rectangle (Fig. 3 B). If the goal is to 
have the computer recognize the shape as a rectangle, pruning 
all branches that appeared as a result of pixelation is essential. 
However, in live-cell imaging, this pruning is not only unnecessary 
but also undesirable. The complexity of the MAT in Fig. 3 D 
accurately captures the complexity of the cell shape in Fig. 3 C, 
where boundary features at all scales are of interest. Thus, the 
high sensitivity of MATs to small-scale features on the bound- 
ary gives our algorithm high accuracy and efficiency for protru- 
sion identification, tracking, and quantification. 

The MAT generates a tree graph that emanates from a sin- 
gle root at the cell center, with each branch ending at a unique 
convex point on the cell boundary (Fig. 3 D). To complete the 
boundary profile, we extend the MAT to the Voronoi graph to 
include concave boundary points as well (Fig. S2 F and Fig. S4). 
With this representation, the distance from the cell center to any 
boundary point can be measured along a path inside the cell. 
This feature turns out to be very powerful for analyzing cell 
shape dynamics because the tips of each protrusion can be easily 
identified as local maxima when the path distance is plotted 
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Figure 2. Boundary artifacts and sensitivity to small boundary perturbations. (A-C] The same cell image with the boundary (cyan) identified at three differ- 
ent threshold values. Bars, 5 urn. (A) Overthresholded case, in which thin filopodia (white arrows) are not properly resolved. Red and yellow lines in A-C 
indicate central lines (skeletons) of protrusions that are longer and shorter than the critical length L cr , respectively. (B) Optimally thresholded case, in which 
all boundary features are properly captured. (C) Underthresholded case, in which filopodia start to merge (white arrows) and a small region of background 
noise (white rectangle) becomes a part of the misrepresented cell boundary. In all cases, FiloTrack identifies and tracks filopodia for a given boundary 
representation, but improperly captured cell features may reduce the accuracy of the final statistics. MovThresh is designed to find optimal thresholds and 
minimize possible boundary artifacts. (D) Two overlaid MATs (black and white) show small perturbation in cell boundaries can cause significant differences 
in the tree graph. Random noise (< 1 -pixel length) creates significant differences in the location of the graph edges and their number so that direct compari- 
son of two MATs is problematic. (E) Rather than comparing MATs directly, we use extracted boundary profiles. This illustrates 1 00 overlaid boundaries that 
were perturbed within 1 -pixel size. Red dots show protrusion tips on each boundary as determined by the corresponding boundary profiles (H). Boundary 
profiles are similar despite differences in tree graphs. (F and G) Expanded views of two regions highlighted in E. The spread of the detected protrusion tips 
(red) is on the same scale as the boundary perturbations, showing that small perturbations do not drastically change protrusion identification when using 
the boundary profiles. (H) 1 00 overlaid boundary profiles of perturbed cell outlines. Red dots show local maxima of the profiles. (I and J) Expanded views 
of two regions highlighted in H and corresponding to the protrusions shown in F and G. Hypersensitivity of the MAT to the boundary perturbation does not 
translate into hypersensitivity of protrusion tip detection using the boundary profiles, even though the latter are created using the MAT. 



as a function of the boundary points (Fig. 3 E). Because pro- 
trusions can have very different geometries (shape, length, area, 
etc.), another advantage of our approach is that it always iden- 
tifies a unique reference point, the protrusion tip. Fig. 4 dem- 
onstrates the accuracy of the Voronoi graph for a cell with 
many protrusions. 

BisectoGraph allows users to import boundaries identified 
in MovThresh or extract them automatically from masks in a 
TIFF stack file. Users can batch process many video files or sin- 
gle frames in a convenient way by specifying file locations and 



numbering in the folder hierarchy. After loading, each boundary 
is regularized to eliminate self-crossings and degeneracies 
(Materials and methods; Fig. S2, D and E; and Fig. S3). Users 
have the option to explore graph structure, resulting boundary 
profiles of the cell, and the length-width relationship of each 
identified protrusion (Fig. S5 A). By interactively specifying 
a critical width for filopodia, users can define the cell body (i.e., 
filopodia-free part of the cell), which can be analyzed in the 
fourth module of the CellGeo, ProActive, which is described 
in the section below on broad protrusions. 
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Figure 3. Using the MAT to identify protru- 
sions. (A) MAT (green) of a rectangle (red). 
The black circles demonstrate that each point 
of MAT is the center of a circle that touches 
the boundary at two or more points. (B, left) 
A tilted rectangle on a discrete grid becomes 
pixelated, creating perturbations along the 
boundary resulting in a complex MAT. Typi- 
cally, to reduce large effects of small bound- 
ary perturbations, pruning algorithms are used 
to remove unnecessary lines, leaving only the 
simple skeleton. (C) Cell with many protrusions 
presents a highly complex shape. Bar, 20 urn. 

(D) MAT of the boundary of the cell in C. In 
this case, the complexity of the MAT prop- 
erly represents the complexity of cell shape. 

(E) Boundary profile of cell in C using the dis- 
tance measured along the paths of the tree 
graph in D. Red dots at local maxima corre- 
spond to tips of cell protrusions. 
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Defining and detecting filopodia-like 
protrusions 

To define, track, and quantify filopodia, we use a very conve- 
nient feature of the MAT: each point of the graph is the center of 
a circle that is tangent to at least two boundary edges (Fig. 4 A). 
The radius of this circle, R, is equal to the shortest distances 
from the graph to the boundary. The tangent circle with the cen- 
ter at the root of the graph and radius 7? max is the largest possible 
circle that can be inscribed in the cell. These properties allow us 
to define geometric scales for cell protrusions. To see this, con- 
sider the distance L from the tip of a protrusion to any point 
along the medial path (the path ending at the graph root). The 
function R(L) is the radius of the inscribed circle at distance L 



and is defined on the interval (0, L path ), m which L pMh is the total 
length of the path. Thus, ^?(0) = 0 and R(L pMh ) = 7? max . 

We can define a critical radius R a < R milK , which determines 
the maximum width of a filopodium (Fig. 4). Using this cutoff 
value, the protrusion length, L pr , is found from the relation 



R, 



(1) 



Using only this criterion, it is possible for the protrusion length 
to be smaller than its width at the base (i.e., L pr < R a ). Thus, to 
ensure only long, thin structures are identified as filopodia, we 
define a minimum length L cl > R CI and require L pr > L cr . Thus, 
a pair of parameters (R a , L CI ) precisely defines a filopodium. 
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Figure 4. Filopodia detection and tracking. (A) Four inscribed circles (a, b, c, and d) with centers on a filopodium skeleton (green). The radius R and dis- 
tance L from tip to center along the skeleton are labeled for the largest circle (d). (A') The shortest distance to the boundary R as a function of the distance 
from the tip t. Vertical lines (a, b, c, and d) correspond to the circles in A. (A") An example of a cell outline (red) and corresponding MAT (green). The 
red box corresponds to the view in A, B, and C. (B and C) Filopodia skeletons (blue, 1, 2, and 3) satisfying a critical radius of 2 or 4 pixels, respectively. 
(B' and C) R[L) curves for medial paths in B and C, respectively. Points where the line for R a intersects with R[L) curves define filopodia bases. (B" and C") 
Resulting filopodia and cell body mapping of the same cell using two different critical radii R cr : 2.5 (B") and 7.5 (C"). (D) An example of protrusion and 
retraction rate analysis for the cell in Video 2. (left) The length of five of the tracked filopodia as a function of time (blue) and corresponding smoothed 
curves (red). Only filopodia that went through the full cycle of protruding and retracting are chosen for analysis, (right) Maximal and minimal velocities of 
filopodia (cyan dots) and corresponding mean values (red lines) calculated using smooth curves in the left graph. (E) The parameter D is a measure of the 
distance between detected filopodium at time t (blue) and f+ 1 (red). If the change in position of two filopodia in consecutive frames is less than D cr , they 
are considered the same filopodium. In this example, tracking only filopodia tips would fail to correctly track filopodia in time. However, FiloTrack correctly 
associates the tips 1 ' with 1 and 2' with 2 (green lines) because D is computed using the skeletons, not just tips. 
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We use the term filopodium skeleton for the fragment of the 
medial path that satisfies these two criteria (Fig. 4, A-C"). 

The critical radius and minimum length are both adjust- 
able parameters in FiloTrack providing built-in flexibility in de- 
fining filopodia. This allows users to customize the definition of 
filopodia to the cell type and analysis to be performed. For ex- 
ample, setting a small critical radius captures only very thin 
filopodia (Fig. 4, B and B'), whereas a large critical radius in- 
corporates wider protrusions (Fig. 4, C and C')- Likewise, ad- 
justing the minimum length defines the threshold for identifying 
filopodia. Scroll bars allow users to quickly and easily set dif- 
ferent combinations of critical radius and minimum length to 
precisely identify filopodia in any cell type (Fig. 5 A). FiloTrack 
automatically corrects for cases when filopodium skeletons 
overlap or when they become wider than the critical width at 
several places along their skeletons (Materials and methods). 

Tracking filopodia over time 

FiloTrack allows users to easily quantify filopodia dynamics, 
measuring length, number, and lifetime. Because FiloTrack 
tracks filopodia over time, it also can compute protrusion and 
retraction rates (Fig. 4 D). Our tracking algorithm is based on 
pairing filopodia from consecutive video frames that are close 
to each other. However, the key feature of our approach is that 
the measure of distance between filopodia is not solely based on 
their tips but on their skeletons (Materials and methods). Thus, 
even when filopodia are too dynamic to make accurate frame to 
frame associations based solely on their tips, FiloTrack still ac- 
curately tracks them (Fig. 4 E). The critical distance, D m defin- 
ing how far a filopodium's position may shift between frames 
and still be identified as the same protrusion, is an adjustable 
parameter in FiloTrack (Fig. 5 B). 

Gap filling. When analyzing videos, filopodia may tem- 
porarily go out of focus or transiently disappear as a result of 
other sources of noise. In these cases, the tracking algorithm 
would record two different filopodia with short lifetimes instead 
of one long-lived filopodia. Well-resolved filopodia might also 
temporarily dip below the critical length L„, and the algorithm 
would treat the reappearing filopodium as a new structure. If 
such events are frequent, statistical measurements will be bi- 
ased toward a high mean filopodia number and short lifetimes. 

FiloTrack was designed to resolve these issues. Once all 
pairs of filopodia from time points t and t + 1 are found, un- 
paired filopodia at time t + 1 are compared with filopodia at 
time t — 1, identifying all filopodia that "blinked" for one time 
step. If unpaired filopodia remain, they are compared with filo- 
podia at time t — 2. The procedure is continued until all time in- 
tervals up to a user-defined time gap, G cr , are tested. Remaining 
unpaired filopodia are considered distinct structures, and the 
algorithm proceeds to the next time frame. The value of the 
maximum time gap G cl can be defined based on biological con- 
siderations. In most cases, G cr set at one to two frames was suf- 
ficient to eliminate the blinking effect without overestimating 
filopodia lifetime. 

Important features of the GUI are that timelines of all 
tracked filopodia are displayed simultaneously on one graph 
(Fig. 5 B, right), and every gap event is indicated as an empty 



circle on the timeline. This provides a quick and easy way to 
assess the quality of the tracking and frequency of gaps. A large 
number of gaps might indicate that either protrusions change 
significantly on a time scale faster than the recording frame rate 
or that the filopodia detection or tracking parameters need to 
be adjusted. 

Lifetime filter. The last parameter built into FiloTrack 
sets the minimum number of frames a filopodium must be pres- 
ent for it to be considered for analysis. The parameter T a may 
be set to 0, and all tracked filopodia will be included. However, 
in many cases, users may choose to only track filopodia with a 
lifetime greater than a specified value. 

FiloTrack allows easy adjustment of the parameters (D m 
G cr , and 7^), so filopodia of any cell type are accurately tracked 
and analyzed (Fig. 4 E, Fig. 5 B, and Video 2). The output of 
FiloTrack includes mean number of filopodia at each time frame, 
number of tracked filopodia, mean filopodial length at each time 
frame, mean length of tracked filopodia, lifetime of tracked filo- 
podia, and SDs from these means. In addition, the time depen- 
dence of the length of all tracked filopodia is saved and can be 
used to calculate the protrusion and retraction rates (Fig. 4 D). 

Quantifying the shape and dynamics of 
broad protrusions 

Lamellipodial size and shape varies greatly between cell types: 
e.g., fish keratocytes have a very large, stable lamellipodium 
encompassing much of the cell, whereas lamellipodia at the 
leading edge of migrating epithelial cells during Drosophila 
dorsal closure are smaller and more dynamic. Additionally, 
cells may produce a mixture of filopodia and lamellipodia, 
using distinct suites of actin regulators. Precise tracking and 
quantification of all types of protrusion dynamics simultane- 
ously in normal, diseased, or genetically manipulated cells is 
required to tease out the protein networks coordinating cell pro- 
trusive behaviors. 

To assess the dynamics of broad protrusions, we again use 
a convenient feature of the Voronoi graph. Every filopodium 
consists of boundary points whose medial paths go through the 
filopodium base. Thus, the remainder of the cell is identified 
simply by excluding these points. This procedure provides a 
computationally straightforward method for partitioning the 
cell into two parts: filopodia and cell body (Fig. 6, A-C). This 
partitioning allows the dynamics of broad protrusions to be 
studied in parallel with filopodial dynamics. 

After removing filopodia, we assess frame to frame changes 
in the remaining cell boundary (Fig. 6, D and D'). The cell 
boundary at time t (Fig. 6, D and D', blue) is overlaid with the 
cell boundary at time t + 1 (Fig. 6, D and D', red). Then, the 
boundary at time t + 1 is compared with that at time t + 2 and so 
on, providing a time-lapse video of cell body dynamics (Videos 
5 and 6). Changes in the boundary are colored white for areas of 
protrusion and black for areas of retraction (Fig. 6, D and D'). 

The ProActive module allows users to display and quan- 
tify protrusion, retraction, or total activity and normalizes re- 
sults to the cell perimeter or total cell area to distinguish 
protrusive activity from changes in total cell area (Fig. 6, E-G')- 
The default in ProActive is to compare successive video frames, 



Computational analysis of cell shape with CellGeo • Tsygankov et al. 449 



^9. HI 




28 

i 26 
£ 

= 24 
1 22 
if 20 
18 



<N> = 


22.48 ±2.25. <L> = 15.04 ± 1.12 ( <L> - L cr 


= 7.44) 










i i 











































10 20 30 40 50 60 70 80 90 




-Filopodia Detection Controls- 



0 Activate Filopodia Detection 
\Jj Show Bisectograph 



critical width = 7.6 



Total Number = 23 
Mean Length = 16.05 
Perimeter = 561 .60 



critical length = 7.6 



- Filopodia Tracking Controls- 



O Activate Filopodia Tracking 



Show Results 

Distance Measure 

Maximum 
t§> Average 
Minimum 



B 




Figure 5. FiloTrack module for detecting and tracking filopodia. (A and B) Screenshots of GUI. The GUI contains a view of the cell, analysis results, and 
control panels for adjusting parameters and viewing results. (A) Filopodia detection. Cell data from BisectoGraph are loaded using the File menu, and the 
cell boundary appears in the cell view section (top left). Checking Show BisectoGraph and Activate Filopodia Detection displays the cell MAT and identified 
filopodia. Red arrows point to scroll bars to select movie frame (top), critical width (middle), and critical length (bottom). Filopodia meeting specified values 
of R cr and t cr are shown in red. Protrusions not meeting these criteria are shown in yellow. The number of protrusions, mean length, and cell perimeter for 
the current frame are highlighted by the red box. The Show Results button loads results for all time frames to the results section (top right) where filopodia 
number (top subsection) and mean length ± SD (bottom subsection) are shown. (B) Filopodia tracking. Clicking Activate Filopodia Tracking and Show 
Results changes the results section to show all time lines of the detected filopodia (vertical green lines). The magenta arrow marks one of the empty circles 
indicating temporal gaps in filopodia tracks. The red arrowhead indicates the scroll bar to select a single filopodium (yellow line in results section and yellow 
dots in cell view [black arrows]). The red box highlights lifetime information for the selected filopodium. Red arrows point to scroll bars for easily changing 
distance, gap filling, and lifetime filter parameters, with cell view and results sections updating accordingly. 
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Figure 6. Detecting and quantifying dynamics of broad protrusion. (A) Cell boundary. Filopodia skeletons are shown in black, and the cell body is shown 
in gray. (B and C) The cell is partitioned into filopodia (B) and cell body (C). (D and D') Changes in cell body shape can be tracked over time. The cell 
body at time t (blue) and t + 1 (red) with regions of protrusion in white and retraction in black. (D') Expanded view of the red box in D. Red arrows high- 
light broad, lamellipodial protrusions. (E) Screenshot of the ProActive module for detecting and quantifying broad protrusions. The cell body with filopodia 
removed is loaded in the left section above a scroll bar for selecting individual frames. The image on the right shows the change in cell body as in D with a 
scroll bar for selecting time lag. Protrusion, retraction, or both can be analyzed and normalized to cell perimeter or area (green curve). A smoothed curve 
through the data is shown in blue. (F) To correct for small cell body fluctuations, the velocity controls (blue arrow) allow users to set a minimum threshold 
for boundary velocity (indicated in red). (G and G') Selecting Run the Range graphs data for the mean and maximum activity for every time lag (G) or for 
every filter value (G'| from 1 to the number typed in the up to box. 
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i.e., with a lag of 1. However, the time lag is an adjustable 
parameter, allowing users to compare cell body outlines at time 
frames that are farther apart in time, e.g., every 5, 10, and 20 
frames. ProActive displays the mean and maximum protrusion 
activity as a function of lag. This provides insight into cell be- 
havior: cells with more persistent protrusions show higher val- 
ues of activity at larger lag values, even if at small lag values 
activity is indistinguishable. Setting the lag to the maximum 
(comparing the first and the last frame) provides a broad view 
of the overall shape change. Finally, if protrusive activity has 
a transient nature, e.g., changes caused by a drug treatment at a 
specific time, ProActive can identify the timing and level of 
peak activity for quantitative analysis of the cells' response to 
the stimulation. 

ProActive also allows users to define the threshold for 
protrusions, allowing small cell boundary fluctuations to be fil- 
tered (Fig. 6 F illustrates exclusion of protrusions moving 8 pixels 
per frame or less). The remaining boundary movements (Fig. 6 F, 
marked in red) represent broad protrusions, such as lamellipo- 
dia. As with other parameters, this filter is fully adjustable for 
accurate detection using any cell type. Identified protrusions 
can then be followed over time and quantified as the fraction of 
the cell edge that protrudes. 

Our approach is also applicable to quantitative analysis of 
cell polarization and motility. The broad protrusion profile 
makes areas of the cell that continually protrude (e.g., the lead- 
ing edge) or retract easily recognizable and quantifiable, making 
our program widely applicable for characterizing both motile and 
nonmotile cells. 

Quantifying growth cone dynamics 

CellGeo also includes the module ConeTrack developed for 
tracking neuronal growth cones. ConeTrack is a natural exten- 
sion of FiloTrack and also relies on the BisectoGraph methodol- 
ogy (Fig. S5 B). The main user-defined parameters in ConeTrack 
that define growth cones are the critical radius of the growth 
cone n ct and neck width W cr > Q„. The critical radius H a is de- 
fined as the minimum allowed radius of the largest circle that 
fits into the growth cone. This parameter sets a minimum size 
for what structures are considered growth cones. The neck 
width W a defines the "base" of the growth cone, that is, the end 
point of the narrow axon shaft between the growth cone and cell 
body. Using these two parameters, ConeTrack automatically 
identifies and tracks every cellular protrusion satisfying the def- 
inition of a growth cone. 

Validation and biological insights 
Filopodial and lamellipodial dynamics. Filopodia are 
thought to sense chemotactic signals, whereas lamellipodia 
drive cell movement. Actin regulators control cell protrusion 
architecture, forming either branched or unbranched actin net- 
works for lamellipodia or filopodia, respectively. Both protru- 
sion types are important for morphogenesis, wound healing, 
and metastasis. 

Many actin-associated proteins have overlapping or re- 
lated functions, and it is unclear what roles they play in regulat- 
ing protrusive behavior in vivo. To address this, we examined 



the roles of Dia and Ena, actin regulators promoting elongation 
of unbranched actin filaments (Goode and Eck, 2007; Bear 
and Gertler, 2009). Dia and Ena each regulate cell protrusions 
in vitro (Gertler et al., 1996; Bear et al., 2002; Peng et al., 2003; 
Gupton et al., 2007; Yang et al., 2007) and in vivo during nor- 
mal development (Schirenbeck et al., 2005; Gates et al., 2007; 
Homem and Peifer, 2008, 2009) and during metastasis (Philippar 
et al., 2008; Sarmiento et al., 2008; Lizarraga et al., 2009). The 
precise roles of Dia and Ena in filopodial behavior in vivo and 
their relationship in this remain controversial. We also investi- 
gated the role of CP, which binds actin filament barbed ends to 
block elongation (Cooper and Sept, 2008). CP can be antago- 
nized by Ena and Dia (Bear et al., 2002; Barzik et al., 2005; 
Schirenbeck et al., 2005; Gates et al., 2009; Hansen and Mullins, 
2010), but the relative contributions of CP antagonism and 
barbed end polymerization by Dia and Ena to protrusion dy- 
namics remains unclear. We thus explored how different combi- 
nations of these three actin regulators regulate the balance 
between filopodia and lamellipodia. 

To do so, we used cultured Drosophila D16C3 cells, 
which naturally form both filopodia and lamellipodia with Dia 
and Ena present in both types of protrusions (Fig. 7, A and B). 
Dia and Ena are the single Drosophila Dia class formin or Ena/ 
Vasodilator-stimulated phosphoprotein family members, mak- 
ing knockdown experiments straightforward. Protrusion dy- 
namics were visualized using fluorescently labeled actin in 
time-lapse videos after overexpression or knockdown of Dia or 
Ena (Fig. 7, C-H'), and the resulting dynamics were analyzed 
using CellGeo (Fig. 8). 

We first elevated levels of active protein by overexpressing 
either a constitutively active Dia mutant lacking the C-terminal 
Dia autoinhibitory domain (DAD; DiaADAD; Fig. 7, D and D') 
or full-length Ena (Fig. 7, E and E') to test the roles of Ena and 
Dia in driving filopodia or lamellipodia. Visual inspection re- 
vealed the expected increase of filopodia as compared with 
controls (Fig. 7 C), consistent with previous work (e.g., Bear 
et al, 2002; Schirenbeck et al., 2005; Gates et al., 2007; Yang 
et al., 2007). CellGeo allowed simple automated quantification 
of filopodia (Fig. 8 A), revealing an increase in both number and 
mean length (Videos 3 and 4). Interestingly, when we used the 
full range of tools in CellGeo, we found a second effect of Ena 
and Dia that was not obvious from visual inspection. Ena over- 
expression (Fig. 7 E') or activating Dia (Fig. 7 D') each signifi- 
cantly increased broad protrusions relative to controls (Fig. 7 C 
and Fig. 8 A). Thus, CellGeo revealed an effect on broad protru- 
siveness that was masked by the dramatic filopodia phenotypes 
(Videos 5 and 6). This highlights the versatility of CellGeo for 
analyzing both filopodia and lamellipodia in the same cell, mak- 
ing it unnecessary to limit analysis to cells with only a single 
protrusion type (e.g., Gupton et al., 2007). 

Next, we used RNAi to knockdown Ena (Fig. 7, H and 
H'), Dia (Fig. 7, G and G'), or CP and examined effects on filo- 
podia and lamellipodia. Although single knockdowns did not 
significantly affect broad protrusions (Fig. 7, F', G', and H'; 
and Fig. 8 B), Ena knockdown reduced filopodia number and 
length (Fig. 7 H and Fig. 8 B) relative to control cells (Fig. 7 F), 
whereas Dia knockdown had no effect (Fig. 8 B). Conversely, 
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Figure 7. CellGeo analysis of genetically manipulated D16C3 cells. (A-A") Drosophila D16C3 cells stained for Dia (A) and phalloidin (A'). (A", inset] 
Dia at filopodia tips (left arrow) and cell cortex (right arrow). (B-B") D16C3 cells stained for Ena (B) and phalloidin (B'|. (B", inset) Ena at filopodia tips 
(right arrow) and cell cortex (left arrow). (C-H) D16C3 cells expressing fluorescent actin and overlaid with the cell boundary (blue) and filopodia (red). 
Protrusions not meeting L cr are in yellow. (C'-H'j Stills. Videos of broad protrusion dynamics (with filopodia removed) for cells in C-H in which every three 
frames are compared (lag 3). Cell bodies at time f (red) or t + 3 (blue) with areas of protrusion in white and retraction in black. (C and C) Wild-type cell 
expressing actin alone. (D and D') DiaADAD expression. (E and E') Ena overexpression. (F and F') Control cell treated with SK RNAi (against BlueScript 
plasmid). (G and G') Dia knockdown. (H and H') Ena knockdown. KD, knockdown; OE, overexpression; WT, wild type. Bars: (C, D, E, H) 5 urn; (A, A" 
[inset], and F) 10 urn; (G) 15 urn. 
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Figure 8. Roles of Ena, Dia, and CP in regulating filopodia and broad protrusions. (A and B) Quantification of filopodia number, length, and broad protru- 
sion activity ( V cr = 4 pixels/frame) of cells in overexpression and knockdown backgrounds. Numbers of cells used for analysis were actin n = 1 6, DiaADAD 
n = 34, Ena overexpression (OE) n = 31, SK knockdown (KD) n = 21, CP knockdown n = 9, Dia knockdown n = 10, Ena knockdown n = 10, Ena knock- 
down DiaADAD n = 7, and Dia knockdown Ena overexpression n = 10. Error bars show 95% confidence interval (mean ± two standard errors). 



CP knockdown significantly increased filopodia number and 
length compared with SK (pBlueScript Sacl-Kpnl) controls 
(Fig. 8 A). Interestingly, CP knockdown also significantly in- 
creased filopodia number and length compared with Ena over- 
expression (Fig. 8 A), suggesting that even with high levels of 
Ena, some actin filament barbed ends remain protected by CP. 
Alternatively, CP knockdown could also affect Dia-based filo- 
podia; however, Dia knockdown did not reduce filopodia num- 
ber or length (Fig. 7 G and Fig. 8 B), suggesting that it is not 
essential for filopodia in these cells. Collectively, our data show 
that Ena is the primary unbranched actin elongator in D16C3 
cells and suggest CP can block barbed end elongation even in 
the presence of high levels of Ena. In addition, Dia is not re- 
quired for filopodia number or length in this cell type but may 
play an alternative role such as regulating cell adhesion (Homem 
and Peifer, 2008). 

Next, we tested the relationship of Ena and Dia in control- 
ling filopodia and broad protrusions by expressing activated Dia 
in Ena knockdown cells or overexpressing Ena in Dia knock- 
down cells. DiaADAD expression in Ena knockdown increased 



filopodial number and length compared with Ena knockdown 
alone (Fig. 8 B). Ena overexpression significantly elevated fi- 
lopodia number in Dia knockdown cells relative to SK controls 
(Fig. 8 B). Interestingly, although individual knockdowns of 
Ena, Dia, and CP did not significantly affect broad protrusions, 
overexpressing Ena significantly increased broad protrusions 
in Dia knockdown cells, and DiaADAD did so in Ena knock- 
down cells (Fig. 8 B). Together, these results suggest that Ena 
and Dia do not require each other to increase filopodia number 
and can act independently in formation of both filopodia and 
broad protrusions. 

Collectively, our data provide new insights into how a sin- 
gle cell type utilizes these three actin regulators. In D 16C3 cells, 
Ena plays the key role in promoting filopodia, with CP antago- 
nizing this activity to balance protrusion number and length. 
Strikingly, CP knockdown is even more effective at elevating fi- 
lopodial number and length than Ena overexpression, suggest- 
ing that blocking CP plays a quantitatively critical role in these 
parameters. Although Dia is not required for filopodia in these 
cells, an activated form can rescue filopodia number and length 
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and promote broad protrusions in cells with reduced Ena levels, 
suggesting that Dia can act independently if activated. We hy- 
pothesize that Ena and Dia are used differentially to create 
structurally similar protrusions with subtle but critical differ- 
ences in dynamics and responsiveness to external cues. Further 
exploration of this hypothesis using CellGeo to analyze how 
Ena and Dia work together to shape protrusions is described in 
detail elsewhere (Bilancia et al., 2014). In summary, CellGeo 
can accurately detect and track filopodia and broad protrusions 
simultaneously in different genetic backgrounds with varying 
morphologies and quantify their dynamics for a precise, unbi- 
ased analysis. 

Neuronal growth cone dynamics 

To demonstrate that CellGeo can be used for different protru- 
sion geometries and different scales of motion and to further 
demonstrate its utility to generate new biological insights, we 
used it to assess the spatiotemporal dynamics of complex pro- 
trusions in neuronal growth cones. Though a large body of re- 
search suggests repulsive cues and collapse factors act through 
the GTPase RhoA (Sebok et al., 1999; Arakawa et al., 2003; 
Yuan et al., 2003; Mueller et al., 2005; Gallo, 2006; Loudon 
et al., 2006; McKerracher and Higuchi, 2006; Woo and Gomez, 
2006; Birkenfeld et al., 2007), little is known about RhoA's role 
in regulating axon extension. Furthermore, growth cones rap- 
idly move and change shape, making it difficult to quantitatively 
assess their dynamic behavior. To overcome this difficulty, we 
developed the module ConeTrack, which allows the immediate 
identification and tracking of growth cones. ConeTrack includes 
user adjustable parameters that can be optimized for each ex- 
periment based on morphological characteristics of the growth 
cones under study. 

To study regulation of growth cone dynamics, we intro- 
duced the RhoA fluorescent activation reporter biosensor (Pertz 
et al., 2006) into catecholaminergic a-differentiated (CAD) 
neuroblastoma cells. Upon serum withdrawal, CAD cells dif- 
ferentiate and extend axons with motile growth cone-like struc- 
tures that phenotypically mimic growth cones of primary neurons 
(Qi et al., 1997). Growth cone collapse and axon retraction can 
be induced by LPA stimulation. We used ConeTrack to analyze 
27 time-lapse videos of CAD cells undergoing axon extension 
followed by collapse after LPA addition (Fig. 9 A). Because 
ConeTrack can simultaneously track multiple structurally dis- 
tinct cellular features, for comparison, we also analyzed cell 
bodies (Fig. 9 B). 

This provided several new insights into Rho regulation 
and growth cone dynamics. The pre-LPA time points revealed 
that active RhoA is already enriched at the edge of both advanc- 
ing growth cones and cell bodies (Fig. 9 C). Strikingly, although 
stimulating the cells with LPA caused the expected global increase 
in RhoA activation, it did not change the shape of the activity 
distribution profiles (Fig. 9 C), with a similar degree of enrich- 
ment at the cortex relative to the interior. In growth cones, peak 
RhoA activation occurred earlier than in cell bodies (Fig. 9 D) 
and on average preceded initiation of growth cone collapse (Fig. 9, 
A and E). After peak RhoA activity, growth cone RhoA activity 
levels rapidly dropped (coincident with collapse), whereas cell 



body RhoA activity was maintained at an elevated level (Fig. 9 D). 
These results provide new spatiotemporal insights into growth 
cone collapse, suggesting elevated RhoA levels in growth cones 
are required to initiate but not to sustain or complete collapse. 
Interestingly, our analysis also revealed a novel, unexpected as- 
pect of growth cone behavior. We found that the time delay be- 
tween LPA treatment and growth cone collapse is inversely 
proportional to the growth cone's size (Fig. 9 F). Although it is 
too early to speculate on the exact mechanisms underlying this 
phenomenon, our analysis provides new insights into axon re- 
pulsion and once again demonstrates the power of CellGeo for 
extracting quantitative information on complex morphodynamic 
behavior not readily obtainable by visual inspection. 

Discussion 

CellGeo provides a highly flexible computational platform for 
defining, measuring, and dynamically tracking cell protrusions 
and other changes in morphology. Our tree graph representation 
of the cell boundary overcomes the computationally challeng- 
ing problem of precisely defining protrusions and provides effi- 
cient methods for quantifying key characteristic features of 
protrusions and cell body dynamics simultaneously. 

We demonstrated the power and versatility of CellGeo by 
using it to analyze the roles of the actin regulators Ena, Dia, and 
CP in cell protrusion dynamics in cultured Drosophila cells, 
revealing their distinct roles in filopodial dynamics. We also 
applied the tool to study growth cone dynamics in CAD neuro- 
blastoma cells and the distribution of RhoA GTPase activity in 
both protruding and retracting growth cones, revealing new in- 
sights into growth cone regulation. 

In its current form, CellGeo is built for applications with 
2D cell images. However, the BisectoGraph methodology is not 
limited to 2D. There are two ways our approach could naturally 
be extended to 3D analysis. The first method involves using a 
z stack to render the cell outline as a 3D surface (or a polyhedron) 
and then building a 3D tree graph inside the resulting 3D shape. 
The second method is to extract 2D features in each z-stack 
plane and "stitch" the extracted 2D features together to form a 
3D representation. We are currently applying the second ap- 
proach to analyze the organization and dynamics of fungal bio- 
films {Candida albicans) in 3D. Imaging in 3D, particularly 
in vivo, will present challenges as a result of the potential decrease 
in signal to noise ratio. This only becomes an issue if it becomes 
too low to allow extraction of the cell boundary. 

Although we focused on protrusions, the computational 
methods in CellGeo are broadly applicable for studying cell 
shape dynamics. For example, we are currently using CellGeo 
to morphologically analyze dendritic spines in neurons, allow- 
ing us to accurately measure and differentiate between the spine 
qualities (thin, stubby, mushroom, and branched) observed dur- 
ing neuronal development. Moving further afield, the method 
can be applied to quantify growth and structural development of 
multicellular formations. For example, we are currently using 
CellGeo to study how the organization and dynamics of fungal 
biofilms influence leukocyte infiltration. Thus, our computa- 
tional method will prove useful to a wide array of scientists 
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Figure 9. Growth cone dynamics and RhoA regulation in CAD neuroblastoma cells. (A) Illustration of growth cone dynamics and RhoA biosensor activ- 
ity in CAD cells as a timeline showing three distinct phases: protruding (top row), reshaping (middle row), and collapse (bottom row). The last two occur 
after LPA treatment. Bars, 10 pm. (B) CAD cell with the growth cone (cyan) and cell body (red) identified by ConeTrack at three time points. Yellow curves 
show the trajectory of growth cone's centroid. Bars, 20 pm. (C) Distribution of RhoA biosensor activity in growth cones (red) and cell bodies (blue) before 
(filled circles) and after (empty circles) LPA treatment. (D) Mean intensity of RhoA biosensor activity at cell edge in growth cones (red) or cell bodies (blue) 
as a function of time. Black dashed line shows the time point when LPA was added. (C and D) Shaded areas show 95% confidence intervals (n = 19 
growth cones; n = 15 cell bodies). (E) Distance from current position of growth cone's centroid to its terminal position, normalized by distance at time of 
LPA treatment (black dashed line). The onset of growth cone retraction (red dashed line) is delayed with respect to time of LPA treatment. (F) Correlation 
between time delay in growth cone retraction and growth cone size immediately preceding LPA treatment (n = 14). Red curve shows a fit to the data using 
the function y = constant/x. 
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studying different biological systems that require automated 
statistical quantification of geometric and dynamic properties 
of the edge. Our user-friendly interface CellGeo allows re- 
searchers to both study geometric aspects of their cell type in 
detail and collect statistical results in an automated and accu- 
rate manner. 

Materials and methods 

Cell culture and immunohistochemistry 

Drosophila D16C3 cells were cultured in Schneider's media supplemented 
with 10% fetal bovine serum, 1% penicillin-streptomycin, and 10 ug/ml 
recombinant human insulin (Gibco/Life Technologies). For overexpression 
experiments, FuGENE HD (Promega) was used to transfect cells with pIZ 
GFP-actin or pIZ monomeric Cherry (mCh)-actin for controls, pIZ GFP-actin 
and pIZ mCh-Ena for Ena overexpression, or pIZ mCh-actin and pIZ GFP- 
DiaADAD for activated Dia expression (mCh constructs were a gift from the 
Rogers laboratory, University of North Carolina-Chapel Hill, Chapel Hill, 
NC). Cells were plated on poly-D-lysine-coated glass-bottom dishes (Mat- 
Tek Corporation) and imaged 48-72 h later. For RNAi, D 1 6C3 cells were 
plated in 1 2-well dishes to 75% confluency and treated with 10 ug of 
double-stranded RNA in fresh growth media on days 1 , 4, and 8. On day 
12, cells were transfected and plated as stated in this paragraph. 

For immunohistochemistry, Drosophila D16C3 cells were plated on 
glass coverslips and fixed in 10% paraformaldehyde in PBS (EM Sciences) 
for 10 min. Cells were permeabilized with 0.1% Triton X-100, blocked 
with 5% normal goat serum, and stained for Ena (1 :500; Developmental 
Studies Hybridoma Bank) or Dia (1:5,000; gift from S. Wasserman, Uni- 
versity of California, San Diego, San Diego, CA). We used Alexa Fluor 
488-conjugated goat anti-mouse or anti-rabbit secondary antibodies 
(1:500; Molecular Probes), and actin was labeled with TRITC-phalloidin 
(1:1 ,000; Sigma-Aldrich). Slides were mounted in Aqua-Poly/Mount (Poly- 
sciences, Inc.). 

CAD cells were cultured in DMEM/F1 2 medium (Gibco) supplemented 
with 8% fetal calf serum, 1% L-glutamine, and 1% penicillin-streptomycin. 
They were imaged in Ham's F-l 2K medium without phenol red (SAFC Biosci- 
ences) supplemented with 15 mM Hepes, 1% L-glutamine, and 1% penicillin- 
streptomycin. To induce differentiation, cells were plated on coverslips 
coated with 18 ug/ml laminin (Sigma-Aldrich). After 90 min, normal cul- 
ture medium was exchanged for imaging medium. Cells were allowed to 
differentiate for an additional 90-1 20 min before imaging. CAD cells 
were transfected 1 2-24 h before imaging with the RhoA biosensor using 
FuGENE 6 (Roche) according to manufacturer's instructions. 3 uM LPA 
(Sigma-Aldrich) was used for LPA stimulation experiments. 

Microscopy 

Imaging of Drosophila D16C3 cells was performed with a lOOx, 1 .4 NA 
oil objective (Nikon) on an inverted spinning-disc confocal imaging system 
(Wallac UltraVIEW; PerkinElmer). Live imaging of GFP- and mCh-tagged 
constructs was performed at room temperature in Schneider's media. 
Images were collected every 2 s for a minimum of 2 min using MetaMorph 
software (Molecular Devices) and a camera (ORCA-ER; Hamamatsu Photo- 
nics). Brightness and contrast were adjusted using ImageJ (National Insti- 
tutes of Health). At least three replicates were analyzed for controls and 
experimental conditions for both live and fixed imaging. Quantification of 
live-imaging videos was performed on 7-34 cells. 

CAD cell imaging was performed on a motorized inverted micro- 
scope (1X81 ; Olympus) equipped with Z-drift compensator focus drift com- 
pensator and total internal reflection fluorescence microscope illuminator, 
a cooled digital 1 2-bit charge-coupled device camera (CoolSNAP; Roper 
Scientific), a 100-W Mercury arc lamp, and MetaMorph imaging soft- 
ware. Images were acquired using a 40x U Plan Fluor N, 1 .3 NA objec- 
tive. A multipass dichroic mirror (89006bs; Chroma Technology Corp.) 
was used with the following band pass filters: CFP (436/20- and 480/ 
40-nm excitation filters) and Forster resonance energy transfer (436/20- and 
540/30-nm excitation filters). Coverslips were mounted in a live-cell cham- 
ber (Attofluor; Invitrogen) and placed on a microscope stage with a heated 
stage adaptor (Warner) that maintained samples at 37°C. Forster reso- 
nance energy transfer/CFP ratio images were calculated after shading 
correction, background subtraction, binary masking, and image registra- 
tion using MetaMorph and MATLAB software (MathWorks, Inc.) as pre- 
viously described (Pertzetal., 2006). Photobleach correction was performed 
with custom MATLAB software that fits the fluorescence decay of the ratio 



intensity over time to a double exponential function and then applies that 
function back to the original images (Hodgson et al., 2006). Time points 
immediately after LPA stimulation were excluded from fitting. 

Boundary representation 

A boundary can be represented as a continuous analytical function 8(x(s), 
y(s)) or a discrete set of connected edges in the form of lines, circular arcs, 
B-splines, conic sections, or other geometric representations. In all cases, 
algorithms exist to construct a medial axis representation of the boundary 
(Fig. 1, C-F). The simplest and most common representation of an object, 
such as a cell, in a binary image is a closed polygonal chain with linear 
segments connecting the centers of neighboring pixels on the boundary of 
the object (Fig. S2, A-C). Such a polygonal chain can be constructed using 
the standard Moore-Neighbor tracing algorithm. In addition to being com- 
putationally efficient and simple, this method preserves the discrete nature 
of the input data and does not involve any preprocessing. However, in 
general, there are two special cases associated with this representation 
that need to be dealt with: self-crossings and degeneracies. 

Self-crossing. Self-crossing refers to the situation in which polygonal 
chain representation of the boundary contains intersections or overlaps 
(Fig. S2, A, D, and E). A simple closed polygonal chain does not contain 
any self-crossings. The MAT of a simple closed polygonal chain is a tree 
graph emanating from a single root with each branch of the graph eventu- 
ally ending at a unique boundary point (Fig. 1, C-F; and Fig. 3 D). We 
refer to a sequence of graph edges from the root to the boundary as a 
path. In contrast, self-crossings create additional features in the resulting 
MAT. For example, a single path might go through several boundary 
points, or some graph edges might overlap with boundary edges. There 
are two options for dealing with self-crossings. One is to keep track of all 
self-crossings and the special features they generate in the MAT at each 
stage of cell shape analysis. The second is to regularize the boundary at 
the beginning of analysis and only deal with MATs in their simplest form. 
We follow the second option. To regularize the boundary, we first elimi- 
nate (set the intensity value to 0) all boundary points with an internal angle 
(i.e., the angle inside a polygon between two adjacent sides) equal to 0 
(Fig. S2 D) and retrace the boundary. This step needs to be repeated until 
no such points are present. Alternatively, one can modify the boundary- 
tracing algorithm, so that unwanted boundary points get eliminated as the 
tracing progresses. Removing boundary points with internal angles of 0 ef- 
fectively removes all boundary extensions of 0 width (i.e., all cell protru- 
sions of 1 -pixel width). All other crossings (such as in Fig. S2 E) can be 
regularized by 1 -pixel dilation of the corresponding boundary pixels (i.e., 
setting the intensity values of their four nearest neighbors to 1). Our solu- 
tion to this problem, which involves adding or removing a single pixel, only 
affects the data locally and on a scale of no more than 1 pixel. 

Degeneracy. Consider for a moment a regular N-sided polygon as a 
boundary. The MAT of this boundary consists of a single vertex at the cen- 
ter and N edges connecting it with the boundary points (Fig. S3 A). The 
single vertex is degenerate in the sense that a small random perturbation 
to the polygon produces a MAT with N — 2 vertices (located near the geo- 
metric center), each with exactly three edges (Fig. S3, B and C). Similarly, 
a simple closed polygonal chain on a regular square grid will likely gener- 
ate degenerate vertices with more than three graph edges. Although, tech- 
nically, this is not a problem: because the MAT is still a tree and the cell 
shape analysis can be performed for the degenerate cases, nondegener- 
ate MATs have well-characterized topological structures and simple geo- 
metric rules for how edges meet at vertices, which we take advantage of in 
our graph-building algorithm (see Tree graph representation). Thus, we 
prefer to deal with nondegenerate MATs after applying random perturba- 
tions to boundary coordinates on a scale of a fraction of the pixel size. 

Tree graph representation 

The tree graph generated by the MAT provides unique paths located inside 
the cell from the root to all convex boundary points but not to the concave 
points (Fig. S2 F). To construct a complete boundary profile with a distance 
measure from the root to all boundary points, the graph needs to be ex- 
tended to include edges connecting concave points to the graph. This can 
be done by considering the Voronoi diagram inside the boundary repre- 
sentation. The Voronoi diagram inside a simple closed polygonal chain is 
a union of the medial axis with all pairs of rays originating from concave 
boundary points at angles perpendicular to the neighboring boundary 
edges and ending on the graph generated by the MAT (Fig. S4). The result- 
ing nondegenerate graph still has all internal vertices of order 3. However, 
strictly speaking, such a graph is not a tree because there are two alterna- 
tive paths connecting the root with each concave boundary point. 
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If the goal is to identify and track cell protrusions, it is sufficient to 
use the tree graph generated by the MAT. However, as we show in the 
previous paragraph, the Voronoi graph is more useful for identifying the 
cell body independent of filopodia. Moreover, our graph construction al- 
gorithm relies on the geometric structure of the whole Voronoi graph, so 
that there is no reason to limit the graph to the subset of vertices defined by 
the MAT. 

Each edge of the medial axis graph bisects the space between two 
boundary components, which are either lines (A, boundary edges) or 
points (p, boundary vertices). Thus, the edges of the MAT consist of three 
types of curves: (1) /n AA , which bisect space between two lines, (2) /i Ap , 
which bisect space between a line and a point, and (3) n w , which bisect 
space between two points. The Voronoi graph also includes edges /m„ that 
connect concave boundary point with the MAT. The curves ,ii AA , /^ pp , and fi,, 
are lines, whereas jii Ap are parabolas. 

When two edges of the Voronoi graph meet/collide at a vertex (an 
event denoted as ®), the third edge emerges (an event denoted as =>) ac- 
cording to the following rules: fn u ® /j, aa => /j, aa (1 ), p AA ® p,. /i v (2), p AA 
® ftp => Map (3), M Af , ® fly => M AA or p, pp (4), /x A( , ® Hk p => Maa or /x pp (5), /j, Ap 
® /x pp => /i Ap (6), p.„ ® (i„ => flip (7), and /i pp ® /u pp => /j, pp (8). The exact 
orientation of the emerging graph edge is defined by the position of the 
boundary components to be bisected, which is determined at the collision 
event, because each colliding edge donates one of two boundary compo- 
nents associated with it to the emerging edge. Thus, all vertices of the Vor- 
onoi diagram can be found in an iterative process starting from the 
boundary and following the collision rules 1-8 until the graph eventually 
converges at the root. Fig. S4 (A-F) illustrates this iterative algorithm and 
the rules 1-7. 

If the boundary consists of N points, k of which are concave and 
(N — k) > 3 are convex, the total number of vertices in the Voronoi graph 
excluding the points on the boundary is N + k — 2. If the boundary points 
are included, the number of vertices becomes 2N + k — 2. When only 
the MAT is of interest, k vertices on the boundary and 2k corresponding 
vertices inside the boundary can be ignored (Fig. S4, D and G) so that 
the remaining MAT consists of N — k — 2 and 2N — 2k — 2 vertices, ex- 
cluding and including boundary points, respectively. Correspondingly, 
the number of edges is 2N + 2k — 3 in the Voronoi graph and 2N — 2k — 3 
in the MAT. Importantly, no matter whether we consider the Voronoi or 
medial axis graph, the tree graph representation requires the storage of 
less than 2N points. To be more specific, it is between N — 2 and 2N — 5 
for the Voronoi graph and between 1 and N — 2 for the MAT, depend- 
ing on the number of concave points, 0 < k < (N — 3). Fig. S4 (G-l) illus- 
trates these numbers with N = 7 and k = 4, so that the MAT (Fig. S4 H) 
consists of four vertices (2N — 2k — 2)— three on the boundary (N — k) 
and one inside [N — k — 2) — and three edges (2N — 2k — 3). The cor- 
responding Voronoi graph (Fig. S4 I) consists of 16 vertices (2N + 
k — 2) — seven on the boundary (N) and nine inside (N + k — 2)— and 
19 edges (2N + 2k - 3). 

Boundary profile 

The tree graph representation of the cell boundary (MAT) has several par- 
ticularly useful features for shape analysis. First, the graph is entirely con- 
tained within the interior of the cell, and each boundary point is connected 
to the tree root by a medial path. This feature provides a mechanism to 
measure distances from an arbitrarily irregular boundary to the cell center. 
Importantly, each local maximum in the plot of such distances, P(n), as a 
function of the boundary point labels (n = 1, 2 . . , N) corresponds to the 
tip of an outward extension on the boundary (i.e., cell protrusion). Protru- 
sions might have very different geometries, but the key step in quantifying 
their geometries is the ability to accurately and robustly identify a single 
reference point, the tip. 

We will be distinguishing two variations of the boundary profiles, 
P(n) and P'(n). The first, P(n), only considers distances measured to the con- 
vex boundary points along the edges of the MAT graph, whereas P'(n) in- 
cludes distances measured from both convex and concave points along the 
edges of the Voronoi graph. Unlike convex points, concave boundary 
points are connected with the root by two different paths so that two values 
of the distance are associated with each concave point. Thus, to deal with 
a single-valued boundary profile, we will use the mean of the two distances 
(alternatively, one can use the smallest of two values). Obviously, protru- 
sion tips should not be sought among concave boundary points. P(n) is suf- 
ficient for that purpose, but P'(n) is a profile of the whole boundary with no 
points taken out of consideration. This becomes useful when, for example, 
one needs to reconstruct the boundary from a truncated graph (see Quanti- 
fying the shape and dynamics of broad protrusions). 



Although, the length of the parabolic segments can be found ex- 
actly, we can simplify and speed up shape profiling by approximating the 
parabolic edges by lines. Formally, this means that we consider exact co- 
ordinates of all graph vertices (as found by a graph-building algorithm) but 
measure distances along the straight lines between the vertices on the cor- 
responding graph paths. This approximation does not alternate the topol- 
ogy of the graph. The shape features that might be affected by this 
approximation are on a scale similar to that of the inaccuracies in bound- 
ary detection and of no interest to us. 

Sensitivity to noise. One of the major disadvantages of the MAT that 
limits its application is its sensitivity to noise. Very small fluctuations of the 
boundary might significantly alter the form of the graph. Fortunately, such 
sensitivity does not affect the important properties of the boundary profile. 
Indeed, a graph perturbed by fluctuations still produces a boundary profile 
with local maxima corresponding to the tips of protrusions. Because we do 
not directly compare graphs, but only the extracted profiles, the sensitivity 
to noise is not an issue (Fig. 2). Again, we note that protrusions with a 
length scale of random boundary fluctuations are not of interest. 

Filopodia identification 

As described in the first section of the Results (see also Fig. 4), a very con- 
venient feature of the medial axis is that by construction each point of the 
graph is the center of a circle that is tangent to two or more boundary 
edges, which serves as a measure of the protrusion width. However, there 
are two easily resolvable complications associated with this approach. 

Uniqueness. As long as P cr < P max , then Eq. 1 in the Results has a 
solution, but it is not necessarily unique. Indeed, a filopodium might be 
slightly wider than the critical width at several places along its skeleton. 
The simplest solution to the nonuniqueness problem is to use for /. pr the 
first intersection of R[L) and R cr . However, in this case, if the typical 
scale for the filopodia width is too close to the critical value R cr , filopo- 
dia might be cut off at length shorter than would be performed by visual 
inspection. Thus, one needs either to choose a bigger value for the pa- 
rameter R cr or to use i pr , which corresponds to the last intersection of R[L) 
and P cr . In the second case, it is useful to introduce a tolerance, e > 0, 
and choose the largest of £ pr 's, satisfying Eq. 1 on the interval from 0 to 
the point where R[L] intersects with P cr + e for the first time. This ensures 
proper identification of filopodia for relatively large cell extensions 
(such as neuron growth cones or spines) that are connected with the 
main cell body by a narrow link. 

Branching. As shown in Fig. 4 C, filopodia skeletons may overlap 
near the base. Thus, if the mean length of filopodia is the quantity of inter- 
est, the contribution from some fragments of cell protrusions will be counted 
twice or more. Once again, this might be an indicator of a poorly chosen 
critical radius P cr . (compare Fig. 4, B and B' with C and C). However, if 
the adjustment of parameters is not an option, a simple way to avoid over- 
counting is to use full-length L p , only for the longest of the partially overlap- 
ping skeletons. For each of the remaining skeletons, the base needs to be 
redefined as the point where the overlap occurs, i.e., at the branching 
point. Because of the tree structure of the MAT, this procedure is well 
defined and straightforward even for the case of consecutively branch- 
ing filopodia. 

Filopodia tracking 

To track filopodia, one needs to specify a distance measure and criteria 
that have to be satisfied to associate a filopodium p at time t with a filopo- 
dium q at time f + 1 . Let rj p = |xP,yPj and ^=(x^,y^J be the skeleton 
coordinates for a pair of filopodia (p, q), in which / =1, 2 . . . N„ and 
(=1,2... N qi and without loss of generality, let the number of skeleton 
points be N p < N q . Then, we define the distance between p and q as 




in which ( )j denotes the mean over index /', and the norm |r| is the Eu- 
clidean distance between two points. Alternatively, instead of the simple 
mean value, one can use the maximum or the minimum value: 

Dp q = max min k p - or 



D = min min r. p - 
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Thus, the matrix D(p,q) tells us how far each filopodium at time f + 1 is po- 
sitioned with respect to each filopodium at time f. 

Introducing a critical distance, D cr , for each filopodium qj(f + 1), we 
find filopodia p m (f) that satisfy the condition 

D( Pm ,q k )<D cr . (2) 

If no such p m is found, filopodium qi is considered to be a newly formed at time 
t + 1 . For the other q/s, we find a pair [p M , q K ] with the smallest distance 

D (PM<%) = minD (P m '%)< 

m,k 

remove the pair from further consideration, and recheck the aforemen- 
tioned condition for the remaining filopodia. This step is repeated until all 
pairs are found. The full procedure is repeated consecutively for each time 
frame so that each filopodium is tracked from the time it appears (i.e., was 
paired for the first time with a filopodium at the next time frame) to the time 
it disappeared (i.e., had no pair among filopodia at the following frame). 
Filopodia that have never been paired are considered to have a lifetime of 
one frame. 

Gap filling. Often there is a need to correct for temporary disappear- 
ance (blinking) of filopodia caused by detection issues. The problem can 
be avoided by filling the gap in the following way. Once all pairs of filopo- 
dia p m (t) and qt(f+ 1) from time points fand f+ 1 are found, unpaired filo- 
podia at time f + 1 are compared with filopodia at time t — 1 on the bases 
of the condition (Eq. 2). This will identify all filopodia that blinked for an in- 
terval of one time step. If there are still unpaired filopodia left, they are 
compared with filopodia at time f — 2. The procedure is continued until all 
gaps up to a desired duration G cr are filled. Only after that, filopodia 
qt(f + 1) that are still unpaired are considered nascent structures, and the 
algorithm proceeds to the next time frame. The value of the maximum time 
gap G cr should be chosen based on biological considerations. If a filopo- 
dium retracts and then a new one is formed in about the same place, it will 
be tracked separately only if G cr is shorter than the time of such reforma- 
tion. Thus, in most cases, taking G cr to be one or two time frames should be 
sufficient to eliminate the blinking effect without overestimating filopodia 
lifetimes. In any case, a pair of parameters (D cr , G cr ) precisely specifies the 
rules by which tracking is performed. 

Statistical analysis 

For data analysis in Fig. 8 and Fig. 9, we used mean values and two stan- 
dard errors (as 95% confidence interval indicated by error bars) of mea- 
surement for individual cells. In Fig. 8, p-values were found using two-sample 
ftest (built-in MATLAB function ttest2). The measurements of filopodia num- 
ber and percentages of protruding cell edge for individual cells were aver- 
aged over the time frames of the time-lapse videos. The measurements of 
filopodia length for individual cells were averaged over all tracked filopo- 
dia during the time-lapse videos. 

Online supplemental material 

Fig. SI shows the MovThresh module for detecting cell boundaries. Fig. S2 
shows examples of boundary representations and self-crossing events. 
Fig. S3 shows examples of degenerate and nondegenerate MATs. Fig. S4 
illustrates the collision algorithm and counting of vertices and edges for 
MATs and Voronoi graphs. Fig. S5 shows the BisectoGraph and ConeTrack 
modules. Video 1 illustrates a MAT interpretation. Video 2 shows an ex- 
ample of a D 1 6C3 cell with a large number of filopodia identified by Cell- 
Geo. Video 3 shows FiloTrack analysis of a control D16C3 cell. Video 4 
shows FiloTrack analysis of a Dl 6C3 cell expressing activated Dia, revealing 
an increase in filopodia number and length compared with control cells. 
Video 5 shows ProActive analysis of the same control D16C3 cell shown 
in Video 3. Video 6 shows ProActive analysis of the same Dl 6C3 cell express- 
ing activated Dia shown in Video 4, revealing an increase in broad protru- 
sions compared with control cells. A ZIP file is also provided that includes 
the CellGeo package, which includes five folders corresponding to the five 
modules: MovThresh, BisectoGraph, FiloTrack, ConeTrack, and ProActive. 
Each folder contains MATLAB codes, data examples, and a text file with 
instructions on how to use the module (input format, output format, and 
processing steps). Online supplemental material is available at http://www 
.jcb.org/cgi/content/full/jcb.201 306067/DC1 . 
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